function lu_decomposition()
% LU_DECOMPOSITION LU分解演示
% 
% 演示LU分解的原理、实现和应用

fprintf('=== LU分解 ===\n\n');

%% 1. 基本LU分解
fprintf('1. 基本LU分解\n');
fprintf('将矩阵A分解为 A = L*U\n');
fprintf('其中L是单位下三角矩阵，U是上三角矩阵\n\n');

% 创建测试矩阵
A = [4, 3, 2; 
     3, 4, 1; 
     2, 1, 3];

fprintf('原始矩阵 A =\n');
disp(A);

% 执行LU分解
[L, U] = lu_decomp_basic(A);

fprintf('L矩阵 (单位下三角) =\n');
disp(L);
fprintf('U矩阵 (上三角) =\n');
disp(U);

% 验证分解
LU_product = L * U;
fprintf('验证 L*U - A 的误差: %e\n', norm(LU_product - A, 'fro'));

%% 2. 使用LU分解求解线性方程组
fprintf('\n2. 使用LU分解求解线性方程组\n');
fprintf('求解 A*x = b，分两步:\n');
fprintf('1) 求解 L*y = b (前向替换)\n');
fprintf('2) 求解 U*x = y (后向替换)\n\n');

b = [9; 8; 6];
fprintf('右端向量 b = [%g; %g; %g]\n', b);

% 第一步：前向替换 L*y = b
y = forward_substitution(L, b);
fprintf('中间解 y = [%.6f; %.6f; %.6f]\n', y);

% 第二步：后向替换 U*x = y
x = backward_substitution(U, y);
fprintf('最终解 x = [%.6f; %.6f; %.6f]\n', x);

% 验证解
residual = A * x - b;
fprintf('残差 ||A*x - b|| = %e\n', norm(residual));

% 与MATLAB内置求解器比较
x_matlab = A \ b;
fprintf('MATLAB解: x = [%.6f; %.6f; %.6f]\n', x_matlab);
fprintf('差异: ||x_LU - x_matlab|| = %e\n', norm(x - x_matlab));

%% 3. LU分解的存在性和唯一性
fprintf('\n3. LU分解的存在性和唯一性\n');
fprintf('LU分解存在当且仅当A的所有主子式非零\n\n');

% 检查主子式
n = size(A, 1);
fprintf('主子式检查:\n');
for k = 1:n
    submatrix = A(1:k, 1:k);
    det_k = det(submatrix);
    fprintf('M_%d = det(A(1:%d, 1:%d)) = %.6f\n', k, k, k, det_k);
end

% 创建一个主子式为零的矩阵
fprintf('\n主子式为零的情况:\n');
A_fail = [1, 2, 3; 
          2, 4, 5; 
          3, 5, 6];
fprintf('矩阵 A =\n');
disp(A_fail);

det_2 = det(A_fail(1:2, 1:2));
fprintf('M_2 = det(A(1:2, 1:2)) = %.6f\n', det_2);

try
    [L_fail, U_fail] = lu_decomp_basic(A_fail);
    fprintf('LU分解成功\n');
catch ME
    fprintf('LU分解失败: %s\n', ME.message);
end

%% 4. LU分解的计算复杂度
fprintf('\n4. LU分解的计算复杂度\n');
fprintf('LU分解的运算量为 O(n^3)\n');

% 测试不同规模的分解时间
sizes = [100, 200, 400, 800];
times_lu = zeros(size(sizes));
times_solve = zeros(size(sizes));

fprintf('\n矩阵规模    LU分解时间    求解时间    总时间\n');
fprintf('--------    ----------    --------    ------\n');

for i = 1:length(sizes)
    n = sizes(i);
    A_test = randn(n, n) + eye(n);  % 确保非奇异
    b_test = randn(n, 1);
    
    % LU分解时间
    tic;
    [L_test, U_test] = lu_decomp_basic(A_test);
    times_lu(i) = toc;
    
    % 求解时间
    tic;
    y_test = forward_substitution(L_test, b_test);
    x_test = backward_substitution(U_test, y_test);
    times_solve(i) = toc;
    
    total_time = times_lu(i) + times_solve(i);
    
    fprintf('%6d      %8.4f s    %6.4f s    %6.4f s\n', ...
            n, times_lu(i), times_solve(i), total_time);
end

%% 5. 多个右端向量的求解
fprintf('\n5. 多个右端向量的求解\n');
fprintf('当有多个右端向量时，LU分解的优势更明显\n');

% 创建多个右端向量
num_rhs = 5;
B = randn(size(A, 1), num_rhs);
fprintf('求解 A*X = B，其中B有%d列\n', num_rhs);

% 方法1：重复使用LU分解
tic;
X_lu = zeros(size(B));
for j = 1:num_rhs
    y_j = forward_substitution(L, B(:, j));
    X_lu(:, j) = backward_substitution(U, y_j);
end
time_lu_method = toc;

% 方法2：直接求解
tic;
X_direct = A \ B;
time_direct = toc;

fprintf('LU分解方法时间: %.6f s\n', time_lu_method);
fprintf('直接求解时间: %.6f s\n', time_direct);
fprintf('解的差异: %e\n', norm(X_lu - X_direct, 'fro'));

%% 6. LU分解的变形
fprintf('\n6. LU分解的变形\n');
fprintf('Doolittle分解: L的对角元为1\n');
fprintf('Crout分解: U的对角元为1\n\n');

% Crout分解示例
[L_crout, U_crout] = lu_crout(A);
fprintf('Crout分解:\n');
fprintf('L矩阵 =\n');
disp(L_crout);
fprintf('U矩阵 (对角元为1) =\n');
disp(U_crout);

% 验证
LU_crout = L_crout * U_crout;
fprintf('验证 L*U - A 的误差: %e\n', norm(LU_crout - A, 'fro'));

end

function [L, U] = lu_decomp_basic(A)
% LU_DECOMP_BASIC 基本LU分解 (Doolittle方法)
% 
% 输入:
%   A - n×n 矩阵
% 输出:
%   L - 单位下三角矩阵
%   U - 上三角矩阵

n = size(A, 1);
L = eye(n);
U = zeros(n, n);

% 检查主子式
for k = 1:n
    submatrix = A(1:k, 1:k);
    if abs(det(submatrix)) < eps
        error('第%d个主子式为零，LU分解不存在', k);
    end
end

% Doolittle分解
for k = 1:n
    % 计算U的第k行
    for j = k:n
        sum_term = 0;
        for s = 1:k-1
            sum_term = sum_term + L(k, s) * U(s, j);
        end
        U(k, j) = A(k, j) - sum_term;
    end
    
    % 计算L的第k列
    for i = k+1:n
        sum_term = 0;
        for s = 1:k-1
            sum_term = sum_term + L(i, s) * U(s, k);
        end
        L(i, k) = (A(i, k) - sum_term) / U(k, k);
    end
end

end

function [L, U] = lu_crout(A)
% LU_CROUT Crout分解
% U的对角元为1的LU分解

n = size(A, 1);
L = zeros(n, n);
U = eye(n);

for k = 1:n
    % 计算L的第k列
    for i = k:n
        sum_term = 0;
        for s = 1:k-1
            sum_term = sum_term + L(i, s) * U(s, k);
        end
        L(i, k) = A(i, k) - sum_term;
    end
    
    % 计算U的第k行
    for j = k+1:n
        sum_term = 0;
        for s = 1:k-1
            sum_term = sum_term + L(k, s) * U(s, j);
        end
        U(k, j) = (A(k, j) - sum_term) / L(k, k);
    end
end

end

function x = forward_substitution(L, b)
% 前向替换
n = size(L, 1);
x = zeros(n, 1);
for i = 1:n
    sum_term = 0;
    for j = 1:i-1
        sum_term = sum_term + L(i, j) * x(j);
    end
    x(i) = (b(i) - sum_term) / L(i, i);
end
end

function x = backward_substitution(U, b)
% 后向替换
n = size(U, 1);
x = zeros(n, 1);
for i = n:-1:1
    sum_term = 0;
    for j = i+1:n
        sum_term = sum_term + U(i, j) * x(j);
    end
    x(i) = (b(i) - sum_term) / U(i, i);
end
end